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Abstract. Variational time discretization schemes are getting of in¬ 
creasing importance for the accurate numerical approximation of transient 
phenomena. The applicability and value of mixed finite element methods 
(MFEM) in space for simulating transport processes have been demon¬ 
strated in a wide class of works. We consider a family of continuous 
Galerkin-Petrov time discretization schemes that is combined with a mixed 
finite element (MFE) approximation of the spatial variables. The existence 
and uniqueness of the semidiscrete approximation and of the fully discrete 
solution are established. For this, the Banach-Necas-Babuska theorem is 
applied in a non-standard way. Error estimates with explicit rates of con¬ 
vergence are proved for the scalar and vector-valued variable. An optimal 
order estimate in space and time is proved by duality techniques for the 
scalar variable. The convergence rates are analyzed and illustrated by nu¬ 
merical experiments, also on stochastically perturbed meshes. 
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1. Introduction 


Numerical simulations of time dependent single and multiphase phase flow and mul¬ 
ticomponent transport processes in complex and porous media with strong hetero¬ 
geneities and anisotropies are desirable in several fields of natural sciences and civil 
engineering as well as in a large number of branches of technology; cf. e.g. [2^ I29j . 
Typically, the discretization in space involves a significant set of complexities and chal¬ 
lenges. MFEM (cf. [T71[^) have proved their potential and capability to approximate 
solutions with high accuracy and physical consistency; cf. e.g. [13 [S]. So far, the tem¬ 
poral approximation of flows and transport phenomena in porous media have received 
relatively little interest (cf., e.g., [23[3[13113031 SllllS] and the references therein) 
and have been limited to traditional non-adaptive first and second order methods, 
even if strong chemical reactions with high temporal variations in profiles are present. 
Rigorous studies of higher order time discretizations are still missing. The low-order 
implicit time discretization is of particular concern with respect to numerical diffusion 
for smooth solutions of transport problems (cf. |45j for a study on numerical diffusion 
for different temporal and spatial discretizations of a transport equation). 

The Galerkin method is a well-recognised approach to solve time dependent problems; 
cf., e.g., [303. However, until now it has rarely been used in practice for discretizing 
the time variable in approximations of initial-boundary value problems. Since recently, 
variational time discretization schemes based on continuous or discontinuous finite el¬ 
ement techniques have been developed to the point that they can be put into use 
(cf. [33 01]) and demonstrate their significant advantages. Higher order methods are 
naturally embedded in these schemes and the uniform variational approach simplifies 
stability and error analyses. Further, goal-oriented error control [3 based on the dual 
weighted residual approach relies on variational space-time formulations and the con¬ 
cepts of adaptive finite element techniques for changing the polynomial degree as well 
as the length of the time intervals become applicable. Variational time discretization 
schemes that are combined with continuous or discontinuous finite element methods 
for the spatial variables are studied for flow and parabolic problems in, for instance, 

[3 H13 0103 03 03 Ell 03 03 03 and for wave problems in, for instance, [30303- 
In these works algebraic formulations of the variational time discretizations are devel¬ 
oped [3 [33133 03133 03, preconditioning techniques for the arising block matrix 
systems are addressed [10313103 and, finally, computational studies are performed. 

Numerical analyses of semidiscretizations in time by variational methods and of vari¬ 
ational space-time approaches can be found in, for instance, [3313103 03 03 . In 
|48j discontinuous variational approximations of the time variable are studied for ab¬ 
stract parabolic problems whereas in |46j their continuous counterparts are analyzed. 
In [33 03 discontinuous variational approximations in time and space are studied and 
error estimates are proved. In m time-dependent domains are considered in an arbi¬ 
trary Lagrangian Eulerian (ALE) framework and the advection-diffusion equation is 
written in mixed form as a system of first order equations in space. In [25] a discontin¬ 
uous Galerkin method in time combined with a stabilized finite element approach in 
space for first order partial differential equations is investigated for static and dynam- 
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ically changing meshes. Error estimates in the L°°{L^) and L^(L^) norm are derived. 
In |34l I35| continuous space-time approximations for nonlinear wave equations with 
mesh modifications and for the Schrodinger equation are considered. Existence and 
uniqueness of the discrete solutions are discussed and error estimates are proved for 
the schemes. 

As far as the MFE approximation of parabolic problems is concerned, in |48| an error 
estimate for the semidiscretization in space is given. However, for the flux variable an 
error estimate is proved only for the norm. No estimate is provided for the error 
in divergence of the flux, that is part of the natural norm of the underlying function 
space Ef(div;H). In [23] and |33| similar error estimates, also in negative norms, are 
presented. In particular, estimates similar to the error estimates for conventional finite 
element approximations are established. The singular behavior of the error estimates 
as t —>■ 0 for initial data in is further included. 

In this work a continuous Galerkin-Petrov (cGP) method is used for the discretization 
in time, whereas the MFEM daEi] is applied for the spatial discretization. Apprecia¬ 
ble advantages of the MFEM are its local mass conservation property and the inherent 
approximation of the flux field as part of the formation itself. In simulating coupled 
flow and transport processes in porous media the flux approximation of the flow prob¬ 
lem is usually of higher practical interest than the approximation of the scalar variable 
itself. To the best of our knowledge, rigorous error estimates for fully discrete varia¬ 
tional space-time discretization schemes that are based on MFE approximations are 
still missing. In our numerical analysis we split the temporal discretization error from 
the spatial one by introducing an auxiliary problem based on the semidiscretization 
in time. We firstly estimate the temporal discretization error and secondly the error 
between the semidiscrete and the fully discrete solution. The order of convergence es¬ 
timates are derived in the natural norms of the variational space-time approach. They 
are summarized in Thm. |4.6| For the scalar variable of the MFE approach one of the 
given error estimates, measured in the norm of L^(0, T; L^(f2)), is optimal in space 
and time if a certain regularity assumption is supposed to be satisfied. For constant 
scalar-valued diffusion coefficients an error estimate for the flux variable in the norm 
of L^(0, T; is further provided. It is optimal in space and suboptimal in time. 

In the Gaussian quadrature points of the temporal discretization optimal order error 
estimates for the flux variable in L^{n) are even obtained for heterogeneous diffusion 
matrices. The existence and uniqueness of the semidiscrete and fully discrete solution 
is further established. Even though a prototype model problem is studied here only, 
we believe that the techniques for analyzing mixed variational space-time approxima¬ 
tion schemes can be applied similarly to more complex flow and transport problems 
in porous media. 

This work is organized as follows. In Sec. our fully discrete variational space-time 
method is developed. In Sec. we address the semidiscrete problem by proving exis¬ 
tence and uniqueness of its solution and error estimates for the semidiscretization in 
time. In Sec. we study the fully discrete problem and show the existence and unique¬ 
ness of its solution. The error between the semidiscrete and fully discrete problem is 
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estimated. In Thm. |4.6| an error estimate for the simultaneous space-time discretiza¬ 
tion is provided by combining the before-given estimates of the temporal and spatial 
discretization. In Sec. we illustrate and validate our derived error estimates by 
numerical experiments. We end our work with some conclusions in Sec. 


2. The fully discrete variational scheme 

2.1. Notation and preliminaries 

Throughout this paper, standard notations are used. Let ft C with d = 2 or 
d = 3, be an polygonal or polyhedral bounded domain. We denote by HP{ft) the 
Sobolev space of functions with derivatives up to order m in L'^{ft) and by (•, •) the 
inner product in L^{ft). Sobolev spaces of vector-valued functions are written in bold 
letters. Further, let = {m £ H^{fl) | m = 0 on dfl} and denote its dual 

space. For the norms of the Sobolev spaces the notation is 

II ■ II := II • llL2(a): II • lip := II • llff'>(a): forpeN,p>l. 

For the mixed problem formulation we use the abbreviations 

V = iT(div; II) = {q £ I V • q £ , W = L^(ft) , 

and 

ll^llv := (ll^f + l|V • 

Let Xq C X C Xi be three reflexive Banach spaces with continuous embeddings. 
Then we consider the following set of Banach space valued function spaces, 

C{I;X) = {u> : [0, T] —)• X I w is continuous} , 

IkWIIx dt < ooj , 

H\l;Xo,Xi) = {wG L2(/;Xo) | dtw £ L‘^iI;Xi)}, 

that are equipped with their naturals norms (cf. [24) 1 and where the time derivative 
dt is understood in the sense of distributions on (0,r). In particular, every function 
in H^{I;Xo,Xi) is continuous on [0,T] with values in X; cf. [21]. For Xq = X = Xi 
we simply write H^{I;X). Moreover, we put Hq{I]X) = {m £ H^{I;X) \ u(0) = 0}. 

For u £ let A : ffo(ft) H~^(ft) be defined uniquely by 

{Au,v) = a{u,v) for all M, z) £ iLg (n) (2.1) 


L^{I-X) = \ w: {0,T) ^ X 


with 


a{u,v) := {DXu,Xv ), 
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where the matrix D = D{x) = {dij{x))f satisfies dij G L°°{^) and is elliptic with 
DMm^>^^D{x)^>D^\^\\ 0Mm^>^^D{x)-^^>9rnm\ ( 2 . 2 ) 


for almost every a; S hi, all ^ € K“ and some constants 0 < Dm < Dm < oo. In (2.2) 


we put 6m '■= Dm ^ and 6m '■= Dm ^■ Under the previous assumptions it holds that 


a(u, u) > a||u ||2 for all u G i?g (fl), 
|a(M,u)| < /3 ||m||i||u||i for all u,v G i?o(^) ■ 


(2.3) 

(2.4) 


Thus, A : Hgifi) >->■ is a linear and continuous operator. For a subspace 

D(A) C i?o(f2) let A : D(A) i—>■ be a bijective linear continuous operator. For 

instance, if fl is a convex polygonal or polyhedral bounded domain and dij G IU^’°°(n), 
for i,j = l,.--d, is satisfied, then the operator A is a bijective linear continuous 
operator from D{A) = n to T^(U); cf. [28] . 

Due to the properties (|2.3|), \2A\ of the bilinear form «(•, •) the lemma of Lax-Milgram 


ensures that the operator A : H^iyi) i—)■ H ^(fl) defined in (2.1) is invertible and 
satisfies in the corresponding operator norm the stability estimates 

ll^ll </? and ||A“^||<a. 

Moreover, for all g G it holds that 

{g,A-^g) = {AA-^g,A-^g) > ^ ■ 

As usual, by c > 0 we denote a generic constant throughout the paper. 


(2.5) 


2.2. Problem formulation 

As a prototype model for more sophisticated multiphase flow and multicomponent 
reactive transport systems in porous media (cf. e.g. [22] |^) we study in this work 


dtu-\7 ■ (DVm) = / 

in n X / , 

(2.6) 

M = 0 

on dfl X I , 

(2.7) 

m(-,0 ) = Mo 

in n , 

(2.8) 


equipped with homogeneous Dirichlet boundary conditions for simplicity only, where 
I = (0, T] with final time T > 0 and the diffusion matrix D satisfies the assumptions 
made in the previous subsection. 


Let / G L'^{I;W) and Uq G Hq{V,) be given. Then the existence of a unique weak 
solution 


u G H^{n)) n Hyi;W) n C{I;W) (2.9) 


to (2.6|-(2.8) is ensured; cf. [221 P- 382, Thm. 5]. We note that (|2.9[ ) already provides 
an improved regularity for the weak solution of (|2.6[)-(2.8); cf. |26l n. 378, Thm. 3]. 
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In order to derive our family of discretization schemes, we first define the auxiliary 


flux variable q := —DVu for the weak solution u of ( |2.6[ )-(2.8 ) that is given by (2.9). 
Since dtu S L'^{r,W) is satisfied by (2.9) and / G L^(/; W) holds by assumption, it 


directly follows that q G L^(/; V). The pair {u, q} G H^{I;W) n C{I;W) x T^(/; V) 
is then also the unique solution to the set of variational equations 

[ {dtu,w)dt+ f {V ■ q,w) dt = [ {f,w)dt, (2-10) 

Jo Jo Jo 

f {D~^q,v)dt— f {u,\7 ■ v) dt = 0 (2-11) 

Jo Jo 

for all w G L'^{I;W) and v G Lp‘{I\V) and satisfies the initial condition m(0) = uq. 


To find (2.11) integration by parts was used. The global problem formulation (2.10) 


(2.11) motivates our semidiscretization in time. 


Remark 2.1. • Below, in order apply Lagrange interpolation in time to the func¬ 

tion f, we need the stronger assumption that f G C([0,r]; W) is satisfied. 

• Below, we introduce a semidiscrete approximation in time of the flux q in a 
subspace of C([0,T];V). For this we need to assume that DVuq G V holds. 


Higher order regularity of weak solutions to (2.6)-(2.8), that is need below for the 


proof of higher order convergence rates, can be obtained under further technical 
assumptions about the data, coefficients and the boundary of the domain LI. For 
the prototype model problem (2.6)-(|2.8|) such higher order regularity results are 
well-known; cf. 


p. 386, Thm. 6]. For (elliptic) regularity results in domains 
with non-smooth boundaries we refer to, e.g., IMW- Below, we tacitly assume 
that the required assumptions about the data and dLl are satisfied such that the 
existence of a sufficiently regular solution can be assumed. Without such an 
assumption the application of higher order methods is not meaningful. 


2.3. Variational discretization in time by a continuous Galerkin 
method 


For the discretization in time we decompose the time interval (0, T] into N subintervals 
In = {tn-ifin]^ where n G {1, ..., N} and 0 = to < < • • • < tn-i < t^ = T. Further 

T denotes the discretization parameter in time and is defined as the maximum time 
step size r = maxi<„<Ar r„, where t„ = — tn-i- We introduce the function spaces 

of piecewise polynomials of order r in time, 


X^{X) 

y^{x) 


G C{I; X) 
G X) 


Ur\J^ G Pr(/n; X) , Vn G {1, . . . , A^}| , 
Wr\I,^ G Pr(/n; X) , Vn G {1, . . ■ ,1V}| , 
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where 


P^(J; X) = \ p-. J ^ X 
and X^{X) C i/i(0,r; W^). We let 


L-.J^X = esX, j=0,...,r| 

3=0 ^ 


X^{X) = ^Ur e A”'(X) w^(0) = 0} . 


Further, we put 


w = Xl{W) X X'^{V) and V = W^(VF) x ^-^(y). 


We equip the function spaces W and V with their natural norms being defined by 

IKWr, <7r}llvv = II^T||i2(/.VK) + ll'5t^T||i2(7.^y) + WItW'l,'^ {I -V) : (2-12) 


i’tIIIv — l|■^T||i2(/.^^A) + ||'i’r||i2(/.v) ■ 

With respect to these norms the space W is a Banach space and the space V is a 
reflexive Banach space. Further, we define the space-time bilinear form Ot € £(>V x 
V; R) by means of 


({^T; Qr }; {^r ; ^r}) 



Wt)) dt 


+ [ {D ^q^,Vr)dt 
Jo 


r 

/ (Ut-, V • Vr) dt 

Jo 


for {ur, Qr} € yj and {wt, v^.} G V. Obviously, the mapping Oi- : W x V M is linear 
and continuous, i.e. 


\ari{Ur,qr}AwT,Vr})\ < c|| || w || {Wr,Vr}\\v 


(2.13) 


with some constant c > 0 independent of r and T. 


For the family of continuous variational time discretization schemes the spaces X^(X) 
of continuous functions act as spaces for the solution whereas the spaces y^~^{X) 
consisting of piecewise polynomials that are discontinuous at the end points of the 
time intervals are used as test spaces. Since the spaces of the trial and test functions 
differ here, a discretization of Galerkin-Petrov type is thus obtained. 


A semidiscrete variational approximation of the mixed form of problem (2.6 1 -( 2.8 1 , ref- 
ered to as the exact form of 00(7-), is then defined by solving the variational equations 
(2.10), (2.11) in discrete subspaces: Find {ur,qr} G X''{W) x X^{V) such that 


/ {dtUr,Wr)dt+ / {X ■ q^,Wr) dt = / {f,Wr)dt, 

Jo Jo Jo 


(2.14) 


{D ^q^,Vr)dt- 


I At ; ^ 
Jo 


' dt = 0 , 


(2.15) 
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for all Wt G ^(W) G with the initial conditions that Wr(0) := Uq 


and q^-iO) := —-DVug (cf. Rem. 2.1). 


We refer to the solution of Eqs. (2.14), (2.15) as the continuous Galerkin-Petrov 


method with piecewise polynomials of order r and use the notation cGP(r). To ensure 
the existence and uniqueness of solutions to (2.14), (2.15), it is sufficient to use the 
test spaces y^~^(W) and y'^~"^{V) with piecewise polynomials of order r — 1, since the 
continuity constraint at the discrete time points tn, n = 0,..., N — 1, that is implied 
by the definition of the solution spaces X'^{W) and yields a further condition. 

By using discontinuous test basis functions u>r(t) = wipn,i(t) and = vtpn,i(t), for 
i = 1,... ,r, with arbitrary time independent functions w G W and v G V, respec¬ 
tively, and piecewise polynomial functions tfn.i : I M. that are of order r — 1 on 
and vanish on /\/„, we can recast the variational equations (2.14), (2.15) as a time 
marching scheme: For n = 1,...,7V find Ut^j^ G Fr{In',W) and G Pr{In',V) 

such that 

J {dtUr,w)ipnA*)dt + j (V • g.^, w) V'n.i(i) dt = j {f,w)ipn,i{t)dt, (2.16) 

J {D~^q.^,v)'tpn,i{t)dt- J (ut-, V • t!) V'n.i(i) dt = 0 (2.17) 

for allw gW and v G V and * = 1,..., r with the continuity constraints Ur | 7 „ (^n-i) = 
(tn-i) and q^| 7 ^(f„_i) = (t„_i) for n > 2 and the initial conditions 

Ur\i„itn-i) ■= Uo, q.,-|7„(in-i) := -DVuq for n = l. 

To determine Ur^j and i represent them in terms of basis functions, with 

respect to the time variable, of the spaces X^{W) and X'^{V) such that 


=J2uiFn,j{t) and q.,^jjt) = J2QiFn,jit), forfeit, (2.18) 
i=o j=o 

with coefficient functions U-f GW and Qi^y for J = 0,..., r and polynomial basis 


functions ipn.j 
points tnj G /, 


r {In 


that are Lagrange functions with respect to r -|- 1 nodal 


satisfying the conditions fn,j{tn,i) = dij for i,j = 0,... ,r. For the 
treatment of the continuity constraint in time we put The other points 

1 ,... r are chosen as the quadrature points of the r-point Gaussian quadrature 
formula on /„ which is exact if the function to be integrated is a polynomial of degree 
less or equal to 2r — 1. The basis functions ipnj G Pr.(/„;K) of (2.18), for j = 
0,..., r, are defined, as usual in the finite element framework, via the affine reference 
transformation onto I = [0,1]. The test basis functions tpn.i G Pr-i{In',F.) with 
'f’n,i{tn,i) = di^i for i,I = l,...,r are defined similarly; cf. [HI |3Z] for details. Now 
we transform all the time integrals in (2.16), (2.17) to the reference interval I. By 


a subsequent application of the r-point Gaussian quadrature formula with weights uJi 
and quadrature nodes ti on I as well as the further notation 


aij := Wi ■ —^f>j{ii) and jdij := (hi ■ 6ij 
at 
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for z = 1,..., r, j = 0,..., r (cf. [HI [351ES])) we obtain the following system of varia¬ 
tional problems for the coefficient functions U^gW and G V of the representation 
(2.18): For n = 1,..., N and j = 1,..., r find coefficient functions {C/^, Qh} G W x V 
such that 


+'r«/ 3 «(V • Qn,w) = Tnfiii{f{tn,i),w) , 
{D-^Ql,v)-{Ul,V ■v) = Q, 


(2.19) 

( 2 . 20 ) 


for i = 1,... ,r and all {rCjt;} G W x V, and where due to continuity in time = 
MT|/„_i(in-i), Q° = qr\i„-i{tn-i) forn>2 and U° := uq, := -DVuq forn = 1. 


Remark 2.2. In the numerical scheme (2.19), (2.20), the flux coefficient functions 


Qi, forj = l,...,r, arise only in the r Gaussian quadrature points tn^i, ■ ■ ■ ■,tn,r G 
{tn-i,tn) of the subinterval Nevertheless, the coefficient functions <5° for n > 1, 
are needed for the unique determination of the semidiserete flux function G X'"(V) 
and an explicit evaluation of q^\j^ by the representation (2.18) in other time points of 
In than in the Gaussian quadrature nodes. The fact that the coefficient functions Qn 


do not arise in (2.19), (2.20) is due to the definition of the Lagrange basis functions 
(Pnj in (2.18) and the fact that the time derivative of the flux variable q does not arise 


in the model equations. 


For the derivation of (2.19), (2.20) from (2.16), (2.17) we tacitly replaced the integrand 
/ on the right-hand side by its Lagrange interpolate Hrf G Fr{In', defined by 


n./(i)i,.=i: f{tn,j)Tn,j{t) iovtGln 

j=0 


( 2 . 21 ) 


We note that the constants fin are satisfying the following property. 

Lemma 2.3 (Coefficient property (C)). There exist constants firm fin G K such that 

0 < fim < fin < fim < oo , for i = l,...,r, (2.22) 

is satisfied. The constants do not depend on the time step size, but only on the number 
r of involved Gaussian quadrature points. 


Proof. Indeed, the coefficients fin = (Li are the Gauss-Legendre quadrature weights 


n 



dx = 


ifi-xlfiP'MW' 


i = 1, 


(2.23) 


scaled to the interval [0,1], i.e. Li = Lij^. 
nomial of degree r and Xi, for i = 1,..., r. 


In (2.23), Pr denotes the Legendre poly- 
are its roots, cf., e.g., [HI p. 436]. Since 
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the sum of the weights Ui equals to two and the weights are all strictly positive, we 
immediately conclude that an upper bound for Cji is given by one. On the other hand, 
we know that \Pl.{x)\ < r{r + l)/2 for any x G [—1,1]; cf. [T^ p. 73]. This gives us the 
lower bound Wi > 2/(r(r + 1))^. ■ 

Below, we will also need the following auxiliary results. 

Lemma 2.4. Let F(t,x) = ^ ^ with coefficient functions 

for j = 0,..., r. Then it holds that 

= I'" (a*F,F)dt=i|lF(t„)f(2.24) 

and 

r 

(2.25) 

i=o 

for some c > 0 independent of Tn- An analogous results holds for coefficients G V. 



Proof. Using the properties of the basis functions ipi and that the r-point Gaussian 
quadrature formula is exact for polynomials of maximum degree 2r — 1 there holds 
that 


{dtF,F)<lt= / / ^ip'^^^{t)Ff^{x)^Lpr,ffit)Ff{x)Ax<lt 


' tn-1 ‘^r2 


i=o 


2=0 


r r pi 


= EE 

j^O 2=0 ' 


dt 


f>,it)-Ut)di{Ff,Fi) 


= E E = E E , K) ■ 

j—0 2=1 2=1 j;=o 


The second of the equalities in (2.241 follows immediately from the first one. It remains 
to prove (2.25). It holds that 


\\F\\lm„;W)<ir+l)J2 / " <c(r + l)5]r„||F;f, 


i=0 ' 


i=0 


with c independent of r^. 


Here we used that ^ ip‘f^ ^{t)dt < CTn] cf. [35l p. 1790]. 


2.4. Discretization in space by the mixed finite element method 

Now, we present the fully discrete approximation scheme that is obtained by discretiz¬ 
ing (2.19), (2.20) with respect to their spatial variables. For this we choose a pair of 
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finite element spaces Wh C W and Vh C V satisfying the inf-sup stability condition; 
cf. [T71 m]. Here, we denote by 7^ = {K} a finite element decomposition of mesh 
size h of the polyhedral domain H into closed subsets if, quadrilaterals in two space 
dimensions and hexahedrals in three space dimensions. Since the software library 
deal, ii [8] that we use for our implementation of the schemes allows only quadrilat¬ 
eral and hexahedral elements, we restrict ourselves to these types of elements in the 
following. Triangular and tetrahedral elements can be treated in an analogous way. 
In our calculations (cf. Sec. we use the Raviart-Thomas element on quadrilateral 
meshes for two space dimensions. For an application in three dimensions based on the 
Raviart-Thomas-Nedelec element we refer to [1511^ . 

The construction of the discrete function spaces Wh and Vh on quadrilateral and 
hexahedral finite elements is done by a transformation Tk '■ K ^ K oi the reference 
element K = [0,1]'^, with d = 2 or d = 3, to the element K through a diffeomorphism 
Tk for all df €Th- We sketch this briefly for d = 2; cf. |1T1[37] for d = 3. For this, let 


:= p: [0,1] 


p(*) = Pi ,> Pi,3 € 

j—0 


We then define the discrete subspaces IF^ C W and C F by 


Wh = WP:=^w€W 


WK°Tf 


-1 

K 


, ^OIK G 


Th} , 


Vh = 


Vl■.= ^^v&v\vKO G gP+i’P X QP’P+I ,iovK&Th]. 


(2.26) 

(2.27) 


The fully discrete continuous Galerkin-Petrov and MFE approximation scheme, re¬ 
ferred to as cGP(r)-MFEM(p), then defines fully discrete solutions Ur,h S X'"{Wh) 
and h ^ X''{Vh) that are represent in terms of basis functions in time by 

r r 

and q^^h\7S^) , for t G, 

i=o i=o 


with coefficient functions t/^ ^ G 


Wh and 


G Vh for j = 0,..., r. The coefficient 
functions are obtained by solving the variational problem (2.19), (2.20) in the discrete 
subspaces Wh C W and Vh C V: For n = 1,... ,N and j = 1,...,r find coeffcient 
functions /j} G Wh x Vh such that 


r 

+r„/3ji(V • Ql,^h^Wh) = Tnfii{f{tn,i),Wh) , (2.28) 

3=0 

{D-^Qlh, Vh) - {Ulh, V-Vh)=0 (2.29) 

for i = 1,... ,r and all {wh,Vh} G Wh x Vh, where G Wh and h ^ h o,re 
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defined by means of the eontinuity constraint (cf. Rem. [S^J, i.e. 


^n.h 'y ' ifn-lj (fn-l) if n > 2 , 

i=o 

r 

Qn,h-= Qn-l.h ^n-l,j{tn-l) if n > 2 , 
i=o 


ul^h ■= PhUo if n = 1, 

Qih-= Phi-DVuo) ifn = l, 


(2.30) 


with Ph : L^{fl) I —Wh and Ph '■ L^{fl) Vh denoting the projections onto Wh 

and Vh, respectively. 


For the derivation of the algebraic formulation of the fully discrete variational prob- 

In 


lem (2.28), (2.29) we also refer to 


the iterative solution of the 


arising linear systems and the construction of an efficient preconditioner is further 
addressed. For solving the algebraic counterpart of Eqs. (2.28), (2.29) we do not ap¬ 
ply an additional hybridization technique as it was done, for instance, in [Hill HD 
and the references therein. We solve the algebraic system by using a Schur comple¬ 
ment technique. In |36j the efhciency of the proposed iterative solver along with an 
adapted preconditioning technique is analyzed numerically. In [TCI [35], the approxi¬ 
mation properties of some families of space-time discretization schemes, including the 
cGP(r)-MFEM(p) approach, in terms of convergence rates and their robustness are 
studied by numerous numerical experiments. Test cases in three space dimensions and 
with heterogeneous and strongly anisotropic material properties are also included. 


3. Existence and uniqueness of the semidiscrete 
approximation and error estimates 


In this subsection we prove the existence and uniqueness of solutions to the semidiscrete 


approximation scheme that is defined by (2.14), (2.15) and its numerically integrated 


counterpart (2.19), (2.20), respectively. The time discretization error is also studied 


in this section. The spatial discretization error is analyzed in Sec. 


3.1. Existence and uniqueness of the semidiscrete approximation 


Theorem 3.1 (Uniqueness of solutions). Let the assumptions of Subsec. 2.2 about 


U, uq and f be satisfied. Then the solution {ut, q^} G X'^iW) x X^{V) of the semidis¬ 


crete problem (2.14), (2.15) is unique. 


Proof. Suppose that {urp.qr.i} G X'^{W) x X^(V) and {uT-,2,Qr.2} G X^(W) x 
X^(V), respectively, satisfy the semidiscrete problem ( |2.14 |, (2.15) and let Ut := 
Ur.i — Ur ,2 and qr := —qr 2 - Then, the tuple {ur^qr} satisfies (2.14), (2.15) with 

/ = 0. Choosing the test function Wr = A~^dtUr G y'^~^{D{A)) C in (2.14) 
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yields that 


f {dtUr,A ^dtUr)dt+ ( (V-q^^A ^dtUr)dt = 0. (3.1) 

Jo Jo 


From estimate (2.5) we get that 

rT 


[ {dtUr,A ^dtUr)dt>-^ f \\dtUr\\H-~UQ)dt. (3.2) 

Jo P Jo 

Using integration by parts in the second of the integrals in (3.1) and recalling that 
A~^dtUr € D{A) C Eq. (3.1) along with (3.2) yields that 

oT nT 

0>^/ \\dtUr\\]j-i(a)dt- I {q^,\/A-^dtUr)dt. (3.3) 

P Jo Jo 


Next, by choosing the test function Vr = DVA ^dtUr G 3^’’ ^(^) in Ep- (2.15) we 
find that 


D q^,DVA ^dtUr)dt 



{DVA ^dtUr))dt = d. 


(3.4) 


Since 


V • {D\JA ^dtUr) = —AA ^dtUr = —dtUr 


and D = by assumption, it follows from ( |3.4[ ) that 

rT 


{Qt^'^A ^dtUr)dt+^ -^||u^fdt = 0. 


Since Ut-(O) = u,-,i(0) — Mr,2(0) = 0 it follows that 


l5tMr)dt+ ^||Mr(r)f =0. 


Combing relations (3.3) and (3.5) shows that 


d>cj liatMrllff-l(a) dt + ^||Mr(T)||^ 


(3.5) 


(3.6) 


and, therefore, Mr = 0. This implies the uniqueness of solutions to (2.14), (2.15). 


To show the uniqueness of solutions q^ of (2.14), (2.15), we choose the test function 
Mr = dtq^ G y''~^{V). Recalling that Mr = 0 by means of the uniqueness result (3.6) 
we obtain from Eq. (2.15) that 


(D ^q^,dtqr)dt = Q. 
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(3.7) 


From (£) ^q^,dtq^) = ^q^,q^) and q^-iO) = 0 we conclude that 


0 =-\\D~^/h,{T)\\\ 


Next, we choose Vr = d^q^ S y'' y'^ C y'' ^(V) by definition, in (2.15), 

recall that rt,- = 0 and use that 


-^^{D ^q^,dtqr) = {D (9tq^) + (D ^q^,d^q^) . 


Together, this implies that 

pT 


0 = J -^{D ^q^,dtqr)dt- J \\D ^^'^dtqrfdt. 


(3.8) 


Since <7^(0) = qr.iid) — <7r.2(0) = 0 and, further, q^iT) = 0 by means of (3.7), it 
follows from Eq. (3.8) along with property (3.7) that q^ = 0. The uniqueness of 
solutions to the variational problem (2.14), (2.15) is thus proved. ■ 


Theorem 3.2 (Existence of solutions). Let the assumptions of Subsec. 2.2 about 
Ll,uo,D and / be satisfied. Then the semidiscrete problem (2.14), (2.15) admits a 
solution {ur,qT.} € X^{W) x A”'(V). 


Proof. To prove existence of solutions to problem (2.14), (2.15), we will use an equiv¬ 
alent conformal formulation, see [43] for a similar approach. 

Find Ur € X’'(77g (fl)) such that Ur{0) = uq and 


/ {dtUr,Wr) dt + / a{Ur,Wr)dt= / {f,Wr)dt 

Jo Jo Jo 


(3.9) 


for all Wr € ^(iJg(fl)). 

The existence and uniqueness of the semidiscrete approximation satisfying (3.9) can 
be established. This is shown in the appendix of this work. Then we define 


Ur '.= Ur and q,. := —D'S/ur. (3.10) 

Obviously, it holds that Ur S X'^(W) since i7o(0) C W. Eurther, we have that 
dfUr G since on each of the subintervals n = 1,... ,N the function 

Ur G X’'(iJQ(0)) admits the representation 

r 

'U'r\I„{t) ='^Ulipn,j{t) , forte In, 

1=0 

with coefficients U-f e 7Jg(0) and polynomial basis functions (pn,j G Pr(/„;K). 
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Next, we prove that S X^iV). 
L'^{r,W) it follows that 


Under the assumption of Subsec. 2.2 that / G 


{-q^,XWr) dt = / (r>Vu^, Vw^) dt = / {f - dtUr.Wr) At / {f,Wr)At 


for all Wt € Y'' ^((^“(U)) with / G L^{I; L^{n)). Thus, we have that 


{-q^,XWr) At = / {f,Wr)At. 


Consequently, it holds that (cf. [Ill p. 18, Eq. (3.38)]) 


[ (S/ ■ q^,Wr)At= [ {f,Wr)At 

Jo Jo 

for all Wt G y’’“^(C'((°(U)) in the sense of distributions. Since / G 

it follows that V ■ qT & T^(/;T^(U)) and, therefore, that G L‘^{I]V) is fulfilled. 

Finally, from the expansion in terms of polynomial basis functions 


qT{t) =-Y^OXW^ipr^At) ^ 

i=o 


(3.11) 


we conclude that q.,. G C'([0,T]; V). 


Eq. (3.91 then directly implies that the functions Ur and q^. defined in (3.10) satisfy 


the first equation of the variational problem (2.14), (2.15). The second equation of 
the system (2.14), (2.15) then follows from the representation ( |3.11 1 of the variable 
qT by testing the identity (3.11) with some function Vt G y^~^{V) and applying the 


divergence theorem of Gauss. Hence, the assertion of the theorem is proved. ■ 

As a corollary of the previous two theorems proving the existence of a unique solution 


to the semidiscrete problem (2.14), (2.15) we obtain an inf-sup stability condition 


within our space-time framework. This result will play a fundamental role in our error 
analyses. For this we need some further notation. Let {uT-,q^} G x X'^{V) 


denote the solution of the semidiscrete problem (2.14), (2.15). We split Ur as 


Urit) = Uq + UtA with U°GA’o(IU). 


(3.12) 


In terms of the tuple {u)!,q^} of unknowns we recast the existence and uniqueness 


result of Thm. 3.1 and |3.2|in the following form. 


Corollary 3.3. Let the assumptions of Subsec. 'L2 about n,uo,D and f be satisfied. 
Let {ur^qr} G X^{W) x X'^(V) be the uniqu e solution of the semidiscrete problem 
(2.14), (2.15) according to Thm. 3.1 and 3.2 Then, the tuple {u°,q^} G A’g(lU) x 


X^{V) with being defined in (3.12) is the unique solution of the following variational 
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problem: Find G Xq{W) x X^{V) such that 


[ {dtU°,Wr)dt+ ( {\/ ■ q^,Wr)dt= f {f,Wr)dt, 

Jo Jo Jo 

[ {D~'^q^,Vr) dt - [ {u°,\7-Vr)dt= f {uo,\7-Vr)dt 
Jo Jo Jo 


for all Wt € y'" ^{W) and v ^(^)- 


(3.13) 

(3.14) 


As a corollary we get the following inf-sup stability condition. 

Corollary 3.4. Let the assumptions of Subsec. \ 2.S\ about Ll,uo,D and f be satisfied. 
Then, there exists a constant 7 > 0 such that 


inf sup 

{«?.q^}GW\{0} {u,^,i,^}gv\{0} 


ari{u°,qr},{wT,VT}) 

IIW,gr}l|w ||{w^,r>r}||v 


> 7 > 0. 


(3.15) 


Proof. The discrete problem (|3.13|) 
Necas-Babuska theorem 


(3.14) satisfies the assumptions of the Banach- 


p. 85]. Since the discrete problem (3.13), (3.14) is well- 
posed according to Corollary |3.3[ the Banach-Necas-Babuska theorem implies the inf- 
sup stability condition (3.15). ■ 


3.2. Estimates for the error between the continuous and the 
semidiscrete solution 


Now we shall show error estimates for the exact form (2.14), (2.15) of the cGP(r) ap¬ 
proach applied to the mixed formulation ( 2 . 10 ), ( 2.111 of our parabolic model problem. 


For this we assume that the following approximation property are satisfied. There 
exist interpolation operators 1^ : Hl{I,W) >-)■ A’q(IF), Jr '■ L‘^{I-,V) X''{V) such 
that for sufficiently smooth functions u G H^{T,W) and q G L{I;V) and all time 
intervals for n = 1 ,..., A^, it holds that 


||m - ) (3.16) 

\\dt{u - IrU)\\L 2 ^^.^w) < Cr;||( 9 [+^M||i 2 ( 7 „.iy) , (3.17) 

Ik - J'r9||L2(/„.v) < cr;+ia[+^q||i2(j^.v) (3.18) 


with some constant c independent of t„ and r. The existence of such approximations 
is obviously ensured, for instance, by using Lagrange interpolation |48) . 

We get the following error estimates in the natural norm of the time discretization. 


Theorem 3.5 (Space-time error estimate for exact form of cGP(r)). Let the assump¬ 
tions of Subsec. 2.2 about LL,uq,D and f be satisfied. Let {u, q} G H^{T, W) xL^{I;V) 
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denote the unique solution of the mixed problem (2.10), ( |2.11 1 that is supposed to be 
sufficiently regular. Then the solution {ur,qr} G x A’’'(V) of the semidiscrete 

problem (2.14), (2.15) satisfies the error estimate 


■-Ur,q- qr}\\w 


1/2 


^ n—1 


where the constant c is independent of Tn, t and T. 


Proof. By splitting 


u{t) = uq + (t) with G Hq{I;W) 


(3.19) 


and recalling the semidiscrete counterpart (3.12), we get that 

u{t) — Ur{t) = u°(t) — u°{t ), dlu{t) = dlu^it) 

for almost every t G (0, T), such that it is sufficient to derive the asserted error bounds 
of the theorem for instead of estimating u — Ur- This will be done in the 

following. 


By (3.16) to (3.18) it holds that 


[u^ - IrU^,q- Jrq}\\w < c|^r2 '-(^||a[+^u°||i2(j^.vv) + 


1/2 


< cT^(^\\dl+\°hm;w) + m+^qhHpv)) ■ (3.20) 

For the discrete functions Wj- := — Itu'^ G d/g (IF), v^. = q.^ — J^q G X^{V) there 

exist, due to the inf-sup stability condition (3.15), functions (pr G Fg (IF), G X'’{V) 
such that 


'y\\{Wr,Vr}\\w\\{TT,i’r}\\v < ar{{Wr, V r}, {pr A r}) 

= ar({u° - IrU°,q- Jrq},WT,fi’r}) 

< C||{U° - IrU°,q- Jrq)\\w\\{TrAr)\\v : 

where the Galerkin orthogonalities 

[ {dt{u°-u°),Wr)dt+ f {\7 ■ {q.^ - q),Wr) dt = 0 , 

Jo Jo 

[ {D~^{q.^-q),Vr)dt- [ {u° - u°,\7 ■ Vr) dt = 0 
Jo Jo 


(3.21) 
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have been used. From (3.21) along with (3.20), we find that 


- IrU^,q^ - Jrq}\\w < n - IrU^,q- Jr9}||w 

< C7“^ (||5[’'’^w||l2(7.w) + \\dl~''^q\\L^(i-v)j ■ 


(3.22) 


From inequality (3.22) along with the interpolation error estimate (3.20) we conclude 
the assertion of the theorem by means of the triangle inequality. ■ 

Thm. |3.5| yields an error estimate with respect to the natural space-time norm of 
the discretization scheme. The estimate is sharp with respect to the contribution of 
\\dt{u — Ur)\\L'^(i;w) to the overall norm (2.12). However, the estimate is suboptimal 
with respect to ||u — Ut\\l‘^(i-w)- In the following theorem, we sharpen our analysis by 
providing an optimal order error estimate also for Ijw — Ur\\L^(i-^w)- This is done by a 
duality argument. For this, the following additional regularity assumption is needed. 

Regularity condition (Rmix)- Suppose that g G L^(I;W). The variational problem, 
find z G H^{T, W) and p G V) with z{T) = 0 such that 


[ {{-dtz,w) + {V ■p,w))dt= [ {g,w)dt, 
JO Jo 

[ {{D~'^p,v) - {z,\/■v))dt = 0 
Jo 


(3.23) 


(3.24) 


for all w G L'^fT, 0; W), v G L^{T, 0; V) admits a unique solution {z,p} G H^{I; L^(H)) 
V) with the improved regularity p G H^{I; V) and the a priori estimate 

\\dtP\\L^{i-,v) < c\\g\\L^(i-^w) ■ (3.25) 

Formally, the corresponding strong form of ( |3.23[ ), ( |3.24 ) is given by 

— dtz + '^-p — g, D~^p + \7z = 0 inHx(0,r), (3.26) 

with z(T) = 0 and homogeneous Dirichlet boundary conditions, that is obtained by 
rewriting the dual problem associated with (2.6)-(2.8), 

-dtZ-V-(Z)Vz) =5 in Hx(0,T), z(T)=0inH, z = 0 on dOx (0, T), (3.27) 

as a system of first order equations. Defining the transformation z{t) := z(T — t) and 
g{f) := z{T — t) we recast (3.27) as a forward parabolic problem in z, 

dtz — V • (DVz) =g in H X /, z(0) = 0 in H , z = 0 on dfl x (0, T), 

such that standard existence and stability estimates can be applied; cf. [26l p. 382, 
Theorem 5]. Then, defining the variable p by means of the second of the identities in 
(3.26), the thus obtained tuple {z,p} satisfies the variational problem (3.23), (3.24). 
Moreover, for g G from [551 P- 382, Theorem 5] we get the a priori estimate 


\\9tz\\L^(I\W) + I1pI1l 2(/;V) < c||g||i2(j.^) . 


(3.28) 
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For this we note that p G L‘^{I] V) can be shown by using the arguments of the proof 
of Thm. 3.2 The a priori estimate of the vector variable p in (3.28) is then a direct 


consequence of the variational equation (3.23). 


Remark 3.6. A regularity condition similar to (Rmiyi) is also used in JM P- 48, Eq. 
(6.16)] to prove the optimal order convergence of a variational time discretization of 
second order parabolic problems in the non-mixed formulation. Currently, it remains 
an open problem how this limiting condition can be avoided in the theoretical analysis. 
The techniques that were developed recently in I25f might be helpful. However, in our 
numerical convergence studies of Sec. the optimal convergence rate that is proved in 
Thm. 3.8 under the condition (Rmix) is nicely observed. 


Below we also need the following auxiliary lemma. 

Lemma 3.7. Let Iq : H^{I, W) ^ y'^{W) and Jq : H^{I, V) y°(y) be interpola¬ 
tion operators that are defined on each subinterval In be means of 

Iou(t) := u(tn-i) and Jov{t) := v{tn-i) for all t G In ■ 

Then it holds that 


\\z — Ioz\\l^{i„,w) < , (3.29) 

\\p — JoP\\l'^{I„-,V) < Tn\\dtp\\L^(I^.V) ■ (3.30) 


Proof. The assertions directly follow from [46l Lemma 6.2]; cf. also [l^ p. 49]. 


Theorem 3.8 {lA Error estimate for the exact form of cGP(r)). Let the assumptions 
of Subsec. \2.1^ about fl,UQ, D and f be satisfied. Further, suppose that the regularity 
condition (Rmix) holds. Let {u,q} G H^{T,W) x L^{I;V) denote the unique solution 
of the mixed problem (2.10), ( 2.11[ ) that is supposed to be sufficiently regular. Then 
the solution {umgr} ^ X''{W) x of the semidiscrete problem (2.14), (2.15) 

satisfies the error estimate 


1/2 


\\u-Ur\\L^(pw) < cr + \\dl+^q\\h^j^,v)) 

U=i 


< CT 


T+1 / 


Proof. We put e„ := G Lfi{I-,W) and Eg := q — q.,. G L'^{I;V) with the 

splitting (3.19) and (3.12) of the scalar variable and its semidiscrete approximation, 


respectively. Further, \ei{z,p} G i7^(0, T; bF)nC'([0, T]; VF) xL^(0, T; V) with z(T) = 
0 denote the unique solution of (3.23|), (3.24) with right-hand side function g = e^. 
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Firstly, recalling that z(T) = 0 and 6^(0) = 0 by definition, we get that 

[ {-dtz.eu) dt =-z{T)eu{T) + z{0)eu{0) + [ (5te„,z)dt= [ {dteu,z)dt. 

Jo Jo Jo 

(3.31) 

Choosing the test function ic = e„ in (3.23) and using (3.31), we find that 

[ ||e„|pdf= f {{dteu,z) + {V ■ p,eu)) dt 
Jo Jo 

= [ {{dteu,z) + {\7 ■ eq,z)) dt- f ((V • e,, z) - (V • p, e„)) df . (3.32) 

Jo Jo 

Choosing the test function v = Eg m (3.24) and recalling that the matrix D is sym¬ 
metric by assumption, we conclude that 


f (S/ ■ eq,z) dt = f {D ^p,eq)dt= ( {D ^eq,p)dt. 
Jo Jo Jo 


(3.33) 


From (3.32) and (3.33) it then follows that 

cT 


[ \\eu\\^dt= f {{dteu,z) + {V-eq,z)) dt- f {{D ^eg,p)-(e„, V-p)) dt. (3.34) 
^0 Jo Jo 


Secondly, by Galerkin orthogonality we find that 

rT 


f {{dteu,Wr) + ■ eq,Wr)) dt = 0 , 

Jo 

[ {{D~^eq,Vr) - ■ Vr)) dt = 0 

Jo 


(3.35) 


(3.36) 


for all Wr € '^r € ^(^)- Choosing Wr = IqZ in (3.35), it follows that 

rT 


[ {{dtCu, Jqz) + (V ■ Sg, Iqz)) dt = 0 . 

^0 


fo 

Further, choosing v^- = Jqp in (3.36) yields that 

rT 


[ {{D ^eq, Jop) - (e„, V • Jop.,.)) dt = 0. 
Jo 


(3.37) 


(3.38) 


Thirdly, combining (3.34) with (3.37) and (3.38), and then using the Cauchy-Schwarz 
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inequality as well as the interpolation error estimates (3.29) and (3.30) yields that 


= / {{dteu, z - Iqz) + {V ■eq,z- Iqz)) dt 

Jo 

- ( ((r>~^eg,p-Jop) - (e„, V • (p-Jop))) dt. 

Jo 

< {\\9teu\\L^{I;W) + \\^q\\L^(I;V))\\z — Ioz\\l'^{I;W) 

+ {\\I^~^h\\eq\\^2(:.v) + l|e«||L2(/;M/))||p - Jop)\\ L 2 {I;V) 

< T{0M\\dteu\\L^(I-W) + \\^q\\L^(I-,V})\\dtz\\L^(i-w) 

+ cr(||eq|li2(7.v) + \\eu\\L^{iiW))\\9tP\\L^(i-,v) ■ 


Applying the a priori estimate (3.28) and the additional regularity assumption (3.25) 


with g = Cu BS well as using the error estimate of Thm. |3.5[ we then find that 

r N 'I 

\\eu\\L^(i-,w) < CT 


< + Wdl+^qWmpv)) ■ 

This proves the assertion of the theorem. 


Next we derive an error estimate for the non-exact form (2.19), (|2.20[) of the cGP(r) 


method. The difference of the non-exact form of cGP(r) to (2.14), ( |2.15 ) comes through 
the numerically integrated right-hand side term in ( |2.19| ). Firstly, we ensure the exis¬ 
tence and uniqueness of the solution to the non-exact form of cGP(r). 


Theorem 3.9 (Existence and uniqueness). Let the assumptions of Subsec. 2.2 about 
LljUo^D and f be satisfied. Then the non-exact form (2.19), (2.20) of the semidiscrete 
problem admits a unique solution {17^, G W xV for j = 1,... ,r and n = 1,... ,N 
defining semidiscrete approximations {urjQr} € X^{W) x A’’"(V) by means of the 


expansions (2.18) and the initial condition 14 ^( 0 ) = Uq. 


Proof. By the definition of the Lagrange interpolation operator 11^ given in (2.21) 
and the representations (2.18) of Ur and q.^ in terms of basis functions we recast the 
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non-exact form (2.191, (2.20) of the semidiscrete problem in the equivalent form 


/ {dtUr,Wr)dt+ / {\I ■ q^,Wr)dt= / {Wrf.Wr) 

Jo Jo Jo 

f {D~^q^,Vr) dt — f (mt, V • Ur) dt = 0 
Jo Jo 


dt 


(3.39) 

(3.40) 


for all Wt S ^(^) and v G ^(^) with the initial condition Mr(0) = Uq. 

Existence and uniqueness of the solution {ur,qr} ^ dl’'(VE) x X''{y) to the system 
(3.39), (3.40) then follows as in Thm. 3.1 and 3.2 with Hr / replacing / in the arguments 
of the proofs. ■ 

Next, we present the corresponding a priori error estimate. 


Theorem 3.10. Let the assumptions of Subsec. \ 2.2\ about LI,uq,D and f be satisfied. 
Suppose that f is sufficiently regular with respect to the time variable. Let {m, q} G 
H^{I;W)xL‘^(I;V) denote the unique solution of the mixedproblem (2.10), (2.11) that 
is supposed to be suffciently regular. Then the solution {wr,9r} G X'^fW) x X^{V) 
of the non-exact semidiscrete problem ( |2.14 1 , (2.15) satisfies the error estimate 

||{u - Ur, q - qAWw <cij2 


+ 'r: 




1/2 


< ct' 


-I- T||9[+^q|||i2(j.vy) -I- r||9[+^/||2,2(/;u/)| • 


where the constants c is independent of Tn, r and T. 

Since the proof of Thm. follows from the proof of Thm. [375| by a standard estimate 
of the interpolation error, we skip it here. For the sake of completeness we summarize 
the proof in the appendix of this work. 


4. Existence and uniqueness of the fully discrete 
approximation and error estimates 


In the first subsection of Sec. we prove the existence and uniqueness of solutions 
to the fully discrete approximation scheme (2.28), (2.29). Then, in Subsec. 4.2 


we 

establish an estimate for the error between the non-exact form of the semidiscrete 


approximation defined by Eqs. (2.19), 


Eqs. (2.28), (2.29). Finally, in Subsec. 4.3 


(2.20) and the fully discrete solution given by 
combine the error estimates of the 


we 


temporal discretization that are derived in Sec. with the error estimates of Subsec. 
|4.2|to get the desired error estimates. 
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4.1. Existence and uniqueness of the fully discrete approximation 


Firstly we prove the existence and uniqueness of solutions to the fully discrete cGP(r)- 


MFEM(p) scheme (2.28), (2.29). For this we need the following lemma (cf. |44l p. 302]. 
Lemma 4.1. For given Wh G Wh there exits a function £Vh satisfying 
\J-Vh=Wh and |l'W/t|| < c||u;/j|| 

for some constant c > 0 depending on 11 and the space dimension d but not on Wh or 
the mesh size h. 

Theorem 4.2 (Existence and uniqueness of solutions). Let the assumptions of Subsec. 
2.2 about Ll,uo,D and f be satisfied. Then the fully discrete problem (2.28), (2.29) 


admits a unique solution {uT,h,QT,h} G X'^iWh) x X''{Vh). 


Proof. Since the fully discrete problem (2.28), (2.29) is finite dimensional and linear 


it is sufficient to show the uniqueness of the solution. The existence is then a direct 
consequence. Assume that there exist two pairs of solutions G X'^{Wh) x 

A’’'(V/i), for fc = 1, 2, that are represented in terms of basis functions by 

r r 

Ur,hit)\i„ = ^n!hTn,j{t) and > for fc = 1, 2 , 


i=0 


i=o 


and t & In with coefficient functions G Wh and & Vh- The continuity con¬ 
straint that is imposed by the definition of X'^fWh) and X^{Vh), respectively, directly 
implies that U°'l = and Further, the pairs ,,(1), q*,,(!)}, for 

k = 1,2, both satisfy the discrete equations (2.28), (2.29). Therefore, it follows that 


- Uf:%Wh)+rnA{y ■ {Q^\ - = 0, (4.1) 

J=0 

(D-HQ^n^ - - {U^nX ” V • ^,) = 0 (4.2) 

for i = l,...,r and all {wh,Vh} G Wh x V/,. Now, by subtracting the equations 
(4.1) and (4.2) from each other and choosing the test functions Wh = U^h ~ 

Vh = TnAiiQhX ~ Qlfh)’ for z = 1,..., r in ([fT]) and (|F^, respectively, we get that 


E -u::i)+rn a,(£»-'(qE- QlXh-QnA = o, (4.3) 

3=0 


for i = 1,... ,r. Summing up Eq. (4.3) from z = 1 to z = r, using Lemma 2.4 and 

rO.l 
^ n.h 


recalling that 1/°’^ = 17°’^ then implies that 


1 


-\\u\Atn) - ulhitn)r + E - QZ)^QhX - Qth) = 0■ (4-4) 
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The symmetric matrix D ^ is positive definite by assumption (2.2) and Pa > 0 under 


the coefficient property (C); cf. Lemma 2.3 Th erefo re, Eq. (4.4) immediately implies 
that for z = 1,..., r. By Lemma 4.1 there exists some Vh G Vh such 


that V ■Vh = — U^t.- Using this Vh as test function in (|4.2|) and noting that the 


ri^'l 


n^h 


'' n^h 


first term in (|4.2|) now vanishes, we obtain that = U^'^, for i = 1,... ,r. This 


n,h'> 


implies the uniqueness of the solution to the fully discrete problem (2.28), (2.29) and 
proves the assertion of the theorem. ■ 


4.2. Estimates for the error between the semidiscrete and the fully 
discrete solution 

In this subsection we derive estimates for the error between the semidiscrete approx¬ 


imation defined by Eqs. (2.14), (2.15) and the fully discrete solution given by Eqs. 
(2.28), ( |2.29 ). Eor this we use the following projection operators (cf. [17], [5] and [40l 
p. 237]) defined in W and V, respectively, by 


Ph-W ^ Wh, {PhW - w, Wh) = 0 


for all Wh € Wh and 
n,. : V 
Ph : V 


Vh, 

Vh 


(V • {UhV - v),Wh) = 0, 

{PhV - v,Vh) = 0 , 


(4.5) 

(4.6) 

(4.7) 


for all Wh S Wh and Vh € Vh, respectively. We point out that Ylh is first defined on 
H^{W) and then extended to V by following [40l p. 237]. For these operators and the 
family of Raviart-Thomas elements on quadrilateral elements for the two-dimensional 
case and the class of Raviart-Thomas-Nedelec elements in three space dimensions 
there holds that 


\\w - Phw\\ < chP'^^\\w\\p+i, 
llt> - n?,^]] < c/i^+^l|t;llp+i, 
||w - Phv\\ < chP~^^\\v\\p+i, 


(4.8) 


||V • (z; - n;,-u)ll < chP+^\\V ■ -ullp+i, (4.9) 
||V • {v - Phv)\\ < c/zP+iJlV • -ullp+i, (4.10) 

for any w £ 77^+^ (fl) and v £ V ■ v £ respectively. 

For the error between the semidiscrete solution and fully discrete we use the notation 

Pu(t) = Ur(t) - Ur,h{t), Eq{t) = q^{t) - Qr^hit) ) 


Pu,n — Eu{tn,i), 


— Eq{^n,i) 


for t £ I and n £ N}, i £ {0,..., r}. Representing the semidiscrete and fully 

discrete solution in terms of basis functions (cf. (|2.18)) there holds that 


Eu(t) = K,n^nAt) and Eq{t) = for i e -f" 


2=0 


2=0 
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Next, we prove two preliminary lemmas. 


Lemma 4.3. Let the assumptions of Subsec. \2.^ about Ll,uo,D and f be satisfied. 
Let the semidiscrete approximation S X^{W) x X^{V) be defined by (2.18)- 

(2.20). Further, let {uT.h j Qr.h } € (Wh) x X''{Vh) be the unique solution of the fully 
discrete problem (2.28), (2.29). Then, for any K = 1,..., N it holds that 


K 


K 


\Eu{tK)f + E E + E E 


n—li—1 


n—1i—1 
K r 

< \\Ur{tK) - PhUr(tKW + C 

71—1 i—1 

(4.11) 

with some constant c > 0 not depending on the discretization parameters h and t. 


Proof. By subtracting (2.28), (2.29) from (2.19), (2.20), respectively, it follows that 

r 

E wh) + r„ A,(V • {Ql - Ql^), wu) = 0, (4.12) 

1=0 

{D-\Ql - Ql^^),vu) - (K - Vh) = 0 (4.13) 

for i = 1,..., r and all {wh, Vh} £ Wh x Vu- For any i = 1,..., r we choose the test 
funct ions w/, = PhU^ - G Wh and Vh = T,,fi.,fiIVhQ\ - Qn,h) G V/t in ( |4.12D 
and (4.13), respectively. By adding the thus obtained equations, using the properties 
of the projection projectors Ph and Ylh defined in (4.5) and (4.6), respectively, and 
summing up from z = 1 to r we get that 


EE PhK - K,h) 

i-l j=0 


(4.14) 


E rn - Qlh), - Qlh) = 0 . 


We note that due to Lemma 2.4 the first term in (4.14) can be rewritten as 

r r 

E E wjiPuUi - PhK - Ulh) 

1=1 1=0 

= 2^\PhEu,n{tn)\\'^ — -\\PhEu^n-l{tn-l)\\'^ ■ 
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Along with some further algebraic manipulations we then conclude from (4.141 that 


2 ~ 2\\PhEu,n-litn-l)\\‘^ 


+ Y,ru - Ql^), Ql - Ql^) 

r 

= Y,rnk{D-\Ql - QIh),QI - n,Q^). 


(4.15) 


i=l 


Recalling assumption (2.2) about D and property (C) (2.22) of the coefficients /Su and 
we obtain from Eq. (4.15) by applying Cauchy-Young’s inequality that 


\\PhEM\\^-\\PhEu{E-l)\\^ +Y,^nPrnBrn\\Qn- 

i^l 
r 

Y,rn\\Qn-nM\ 


< 


Ph 

Pm 


(4.16) 


i=l 


Summing up inequality (4.16) from n 
shows that 


1 to AT and noting that PhEu{to) = 0 then 


\\PhEMf + EErnPM-Ql^l 

n—1 i—1 


< 


P 


K r 


(4.17) 


'm 




n=l2=1 


for any AT S N with K < N. By using now Lemma |4.1[ there exists for any 
i G a. Vh e Vh such that V ■ Vh = PhE^,^ and |lt;,i|| < c\\PhEl^J\. By 

testing (4.13) with this v^, we get by using the Cauchy-Schwarz inequality along with 
assumption (2.2) about D that 

\\PhElJ<ceM\\Ql,-Qln\\, (4.18) 

for n = 1,..., iV, i = 1,..., r. 

Combining (4.17) with (4.18) is follows that 


i\ T 

\\PMtK)r + E T.^nPm\\Q\ - Q\,h 

n—1i—1 

K r 


K r 

n—1i—1 


(4.19) 


,2 ||2 


12 = 1 2=1 
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By = Ql- Q'l ^ and the triangle inequality relation ( |4.19[ ) implies that 

K r K r 

\\Eu{tKW + E E + E E ^nWElS 

n—1 i—1 n—1 i—1 

K r 

<cEE-^ IIQ^ - UhQnW'^ + c||Ur(tfc) - PhUritkW 

n^l i-1 
K r 


(4.20) 


n—1i—1 


Observing that E]^^ — PhEl^ ^ = JJ^ — PhU^, inequality (4.201 proves (4.11). 


In the second lemma we restrict ourselves to the case that D = dl with some d > 0 is 
satisfied. An extension of the provided estimates to more general matrices D{x) still 
remains an open problem. 


Lemma 4.4. Let the assumptions of Subsec. \2.S\ about Ll, uq and f be satisfied and D = 
dl with some d > 0. Let the semidiscrete approximation {ur,qr\ G X^^iW) x X'^{V) 
be defined by (2.18)-(2.20). Further, let {ur,h, QT.h} ^ X^{Wh)xX'"(Vh) be the unique 
solution of the fully discrete problem .28[ ), ( |2.29[ ). Then, for any K = . ,N it 

holds that 


K 


n—1 i=l 


E E + WPhESkW 

K r 


(4.21) 


n—1 i—1 


Proof. Introducing the projectors into the error equations (4.12)-(4.131 yields that 

(4.22) 


^ ^ ^ij {PfiEilj^ j,^, Wfi) Tn fin (V • yj, Wjfi — 0 , 

1=0 

{PuEl^^,Vk) - {PhEl^,V -Vk) = Q 


(4.23) 


for n = i = 1,..., r and all {wh, Vh} G Wh x Vh- Observing that for any 

n > 2 the quantities and „ are linear combinations of and 

for i = 0,... ,r, respectively, and that PhE^ i = 0 and PhE^i = 0 by definition of 


{Ui,Qi} and {U^ i,,Qnh}, it follows that Eq. (4.23) is also satisfied for j = 0 and 
any n > 1. Using this, we obtain by multiplying (4.23) with ctji and summing up the 
resulting identity from i = 0 to r that 


\j=o / \i=o / 


(4.24) 
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for any & Vh- We note that we changed the notation for the indices. By testing now 
( |4.22[ ) with Wh = Ej=o ^ij^hEi n e Wh and ( |4.24[ ) with Vh = G V?,, we 

get by summing the resulting equations and using the inequalities of Cauchy-Schwarz 
and Cauchy-Young that 


I _ 2 ■^ / 

'y ' ^ijEhE^ ,^ + Tn Pii( ^ ^ ^ijPhEq^rn EhEq^ 


f=0 


J=0 


= T„ Pul ^ aijPhEi,^, V • (P/i - 


j=o 


1II 

— 9 ^^EhEi ,^ 


j=o 


\lr^PU^-{Ph-nn)ElJ 


for n = 1,..., iV and i = 1,..., r. The inequality above further simplihes to 

I ^ijEhEl ,^ + 2r„ /3ii/ '^Q.jjPhEii ^^ Pj^E^ ,^ 

0=0 0=0 (4.25) 


<t^PI\\W-{P^-U^)E 


q,n\\ ’ 


for n = l,...,iV, i = l,...,r. Dividing (4.25) by Pu (note that pu > 0 for all 
i = 1,..., r), summing up the resulting inequality from z = 1,..., r and using Lemma 
|2.4| gives that 


o II ^ijEhEl^ ^ + \\PhEq(tn)\\ <\\PhEq{tn-l) 

i=l '''n Pu j=o 

r 

+ ■{Ph-n^)ElJ 


(4.26) 


forn = 1,..., A^. By summing up (4.26) from n = 1,..., K and noting that PhEq(tQ) = 
PhE^ 1 = 0 for the choices of the semidiscrete and fully discrete coefficient functions 
and Q^ji (cf. their definition below (2.19), (2.20) and Eq. (2.30l) we get that 

K r 


EE-Vl^ aijPhEl^^ 

n=li=l '’’n Pu 


0=0 


\PhEq{tK)\Y 


K 


(4.27) 




n—1 i—1 


We now estimate the divergence of the flux. By testing (4.221 with Wh =y -II/iElq „ G 
Wh, and using the inequalities of Cauchy-Schwarz and Cauchy-Young {pu > 0 for all 
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I = 1,..., r) we get that 


IV • = -( E ^^jPhEi^n, V • n^El^ 

j=0 

1 II T B 

-L II ^ ^ , 'n H', 


< 


^Ttt, I3ii 




j=0 


V • UhE 


q,n 


for n = 1^..., N, i = 1,... ,r. Summing up the previous inequality from n = ..., K 

as well as from i = 1,..., r, using (4.271 along with {Ph — 'n.h)E^^^^ = (Ph — Tlh)Qh 


.niP, 


by definition of the projectors Ph and II/j we obtain that 

K r K r 

E E • ^hE\J^ < E E V • {Pn - Tlu)Q\ 

n—l i—1 n—1 i—1 

which proves the assertions of the lemma. ■ 

Now we combine the inequalities of the previous lemmas to estimate the error between 
the semidiscrete and the fully discrete solutions in the norms of LB {I ]W) and LB{I\V). 


Theorem 4.5. Let the assumptions of Subsec. \2.S\ about n,uo,D and / be satisfied. 
Let the sufficienctly regular semidiscrete approximation {ur, q^} S X^{W) x X^lfV) he 
defined by (2.18)-(2.20|. Further, let {ur,h,qT,h} ^ x X'"{Vh) be the solution 

of the fully discrete problem (2.28), ( |2.29[ ). For the scalar variable Ur it holds that 

\\ur — UT^h\\L^(I-,W) ^ ch^^^. (4.28) 

For the vectorial variable q.,. it holds that 

/ N r 


. 1/2 

2 1 < chP+^ . 


( E E - qrAtn,i 

\n^l / 

Further^ for D{x) = dl, for some d > 0, it holds that 

WQt - QT,h\\L^i-,L^(n)) < chP+^ 

and 

/Nr \ 

( E E - qr.h(.tn,i)\\v ) < chP+^ . 

\n —1 2—1 / 

The constant c does not depend on the discretization parameters h and r. 


(4.29) 


(4.30) 


(4.31) 


Proof. By using (2.25) and recalling that = Euitn-i) we find that 

.Nr N . 

\\nr - ^ d E E Pr^\\ElJ^ + Yi r„||K(t„-i)f j . (4.32) 

^ n—l i—1 n—l ' 
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By inequality (4.111 we get for the first term on the right-hand side of (4.32 ) that 


N 


N 


E E II II '+ E E II f <c(J2J2'^JWn-PhK 


N 


n=l i=l 

N r 


i\ r \ 

+ TN\\Ur{tN) - PhUritNW + E ) ■ 

n=l i=l ^ 


(4.33) 


Using (4.111 again, we find for the second term on the right-hand side of (4.32) that 

\\Eu{tK)\\'^ < \\Ur(tK) — PhUT{tk)\\^ 

^ ^ (4 34) 

+ ^E Edl*^-- ^^'^"ll'+ ll<3n - ) 

n— 1 2=1 


for iiT = 1,..., A^. Combining now (4.32) with (4.33) and (4.34) and using the approx¬ 
imation properties (4.8)-(4.10) of the projection operators we then get that 


— Ur,h\\L^{I-W) 


N 




2 = 0 


X+i 

N 


(4.35) 


-b ^max^ ||u^(tK)||p+i + E Ell Qnllp-l-l j ) 

n —1 2=0 


where the arising constant does not depend on the discretization parameters h and r. 
The result (4.28) directly follows from (4.35) under the assumption of the theorem of 
sufficiently regular coefficient functions € W x V. From (4.33) along with 

(14.341) and (|4.8|)-(|4.10|) we further conclude that 


N 


E E ll^q(^"4)f < . 


(4.36) 


n —1 2=1 


This proves (4.29). 


By using (2.25), recalling that = Eq{tn-i) 
projection operator we hnd that 


and applying the boundedness of the 


WQt 9r.?illL2(7.£,2(f2)) 

< c(\\Qt ~ ■P?i9rllL2(/;i2(n)) + WPhQr ~ 9t.?i II L 2 ( 7 . 7,2 (q)) ^ 

.Nr N 

^ d E E ^nWElX + E rr.\\PhE^{tu-iX 

^n— 1 2—1 n —1 

+ WQt - -P/i9rllL2(/:I,2(Q)) ) • 


(4.37) 
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For D{x) = dl the second term on the right-hand side of (4.37) can be bounded 
from above by means of the inequality (4.21) along with the observation that {Ph — 
„ = {Ph — 'i^h)Q^n by definition of the projectors Ph and Il/i. Recalling further 
the boundedness of jda (cf. Lem. 2.3) we conclude that 


K 


\PhE^{tK)r < cY,Y.^n\\{Ph - n,)Q^ 


(4.38) 


n—1i—1 


for iF = 1,...,A^. Finally, combining (4.37) with (4.33) and (4.38) and using the 
approximation properties (4.8)-(4.10) of the projection operators we then get that 

Ikr 


N 


< C 


||[/*||p+i ^max 

n—1 


2=0 


N 




(4.39) 


N 


n—1 2=0 


where the arising constant does not depend on the discretization parameters h and 
T. The result (4.30) directly follows from (4.39) under the assumption of sufficiently 
regular coefficient functions {Ul,Ql^} gWxV. 


To estimate the divergence part of the error in (4.31), we use that by definition of the 
projection operators it holds that 


N 


11^ ■ 


n—1 i—1 

N 


N 


(4.40) 


< E E11^ • + E E11^ • 


n—1 


i-1 


n—1 


2=1 


The assertion (4.31) then follows from (4.40) combined with (4.21) and the approxi¬ 
mation properties (4.8)-(4.10). ■ 


We remark that the inequalities (4.30) and (4.31) provide an error control for the 
spatial discretization in the Gaussian quadrature points or temporal degrees of freedom 
of the subintervals In with respect to the norm of T^(n) and V, respectively. For an 
error control with respect to the norm oi {I; V) or L'^{I-,V) a further estimate of 
ifq „ is required which remains an open problem. 


4.3. Error estimates for the error between the continuous and the 
fully discrete solution 

In this section we combine the results of Thm. 13.81 and Thm. 13.101 with the estimates 
of Thm. [T5]to prove the convergence of the fully discrete scheme. 
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Theorem 4.6. Let the assumptions of Subsec. \2.S\ about LI,uq,D and / be satisfied. 
Let {u,q} S H^{L;W) x L^(J; V) denote the unique solution of (2.10), (2.11) that is 
supposed to be sufficiently regular. Further, let {uT,h^ i T.h} € T''(V14,) xX'^{Vh) be the 
uniquely defined solution of the fully discrete problem ( |2.28[ ), (2.29), respectively. Sup¬ 
pose that the semidiscrete problem (2.19), (2.20) admits a sufficiently regular solution 
{ungT-} G X^fW) X X^{V). Then, there holds that 

||w - Ur,h\\L^I-,L2{Q)) < c('r’' + . (4.41) 

For homogeneous diffusion coefficients D{x) = dl, with some constant d > 0, there 
holds that 

Wq - < c{t^ + hP+^). (4.42) 

Under the regularity condition (Ry^m) given in ( 3.25[ ) and for interpolated right-hand 
side functions ( 2 . 21 ) there holds that 


' - Mr,/i||L2(/;i2(a)) < c(t’'+^ + /lP+^) . 


(4.43) 


The constant c in (4.41 )-( 4.43[ ), respectively, does not depend on the discretization 
parameters h and r. 


Proof. By using the triangle inequality, Thm. [3T0| and Thm. [4)5| it follows that 
||u — < 2\\u — MT||i 2 (/.i 2 (Q)) + 2\\Ur — 11 ^2 ( 7.^2 (Q)) 

<c(r2'- + /72(P+i)) , 

where sufficient regularity of the continuous and semidiscrete solution with appropriate 


upper bounds for the solutions (cf. Thm. 3.10 and Thm. 4.5) is assumed. The inequality 
(4.42) is obtained similarly. The estimate (4.43) can be concluded in the same way by 


using now the result of Thm. |T8| instead of Thm. |3.10| 


Remark 4.7. 


The error estimate (4.43) is optimal in time and space. The 


assumption of an interpolated right-hand side function ( 2 . 21 ) can still be dropped 


even though this is not explicitly done in this work. It requires to estimate the 


error between the exact form of cGP(r) defined in (2.14), (2.15) and the fully 
discrete solution. In this case the arguments used to prove Thm. \4.I\ have to 
be augmented by an estimate of the interpolation error for the right-hand side 
function, similarly to the proof of Thm. \3.1C\ 


The error estimate (4.42) is suboptimal in time. It remains an open problem 


to analyze if the estimates can still be sharpened to order r + 1. In our nu¬ 
merical study presented in Sec. convergence of order r + 1 will be observed 
for the temporal discretization of the scalar and the flux variable. Moreover, 
this is even observed in the (spatially) stronger norm oflf{Q,T\V) instead of 
Lfff),T-,Lfl{VL)) for the flux variable. 
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5. Numerical studies 


In this section we present numerical studies in order to illustrate the error estimate 


given in Thm. 4.6 for the fully discrete scheme (2.281, (2.29) combining a variational 


time discretization with the MFEM. Moreover, we analyze the robustness of the con¬ 
vergence behaviour with respect to random perturbations of the meshes. Thereby 
we mimic mesh distributions of applications that are of practical interest. Additional 
convergence studies for variational space-time discretizations of the proposed type as 
well as for discontinuous time discretizations can be found in [UlISZl for parabolic 
problems and in [saisT] for variational space-time discretizations of wave equations. 
In [37l [15] the efhcient iterative solution of the resulting algebraic system of equations 
(2.28), (2.29) along with the construction of appropriate preconditioning techniques 


is carefully addressed. In the literature, further computational studies of variational 
time discretization schemes are presented also for different kind of flow and transport 

problems in, e.g., n UllSlIli I3H1I311301131] • 

In order to determine the space-time convergence behavior we consider in our numerical 
study the cGP(2)-MFEM(2) approach. That is ( 12.14 |-(2.15) with r = 2 combined 
with the mixed finite element method MFEM(2) based on the choice p = 2 in the 
definition ([2.26 1 and (2.27) of the tuple of MFE spaces. We prescribe the solution 


u^ix, t) := sin(a;f) sin(7ra;i) sin(7rx2), in x (0, T), 


with n = (0,1)^, u! = IOtt of problem (2.6)-(2.8). The corresponding flux function 
is then given by = —DWue for D — I. We choose the final time T = 1. On 
the coarsest level (level 0) the temporal mesh is uniformly refined into = 10 time 
subintervals and the corresponding spatial mesh consists of a single cell. To determine 
the experimental orders of convergence the space-time mesh is refined uniformly by a 
factor of two in each of the space dimensions and in the time dimension. We use the 
abbreviation 


■= UEit) - Ur,hit) and := q^^it) - qr,hit), 

where we denote by Ur,h and by q,. ^ the fully discrete cGP(2)-MFEM(2) approxima¬ 
tion of the primal variable and the flux field. The discretization errors for are 

measured in the L^(/; L^(n))-norm and for in the L^(/; V)-norm. As usual, 

the integral over the spatial domain 17 and the integral over the time domain I = (0, T) 
in the error norms are evaluated elementwise in space and time by appropriate quadra¬ 
ture rules of sufficiently high order of accuracy. 


5.1. Uniform meshes 


We summarize the calculated errors and their experimental order of convergence 


(EOC) for the proposed space-time discretization in Tab. 5.2 and further illustrate 
them in Fig.|5.1[ The numerical results confirm the expected third order rate of conver¬ 


gence established in Thm. 4.6 (cf. also Rem. 4.7) for the discretization in the space-time 
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Level 

N 

Tn 

\%\ 

h 

A^DoF 

0 

10 

l.OOOe-01 

1 

1.4142e-00 

33 

1 

20 

5.000e-02 

4 

7.0711e-01 

120 

2 

40 

2.500e-02 

16 

3.5355e-01 

456 

3 

80 

1.250e-02 

64 

1.7678e-01 

1776 

4 

160 

6.250e-03 

256 

8.8388e-02 

7008 

5 

320 

3.125e-03 

1024 

4.4194e-02 

27840 


Table 5.1: Space-time mesh with number of time subintervals N, global time discretiza¬ 
tion parameter r„, number of cells \Th\, global space discretization param¬ 
eter h and degrees of freedom A^doF per degree of freedom in time. 


Level 

II cGP( 2 )|| 

lr“ llL 2 (/;L 2 (n)) 

EOC 

II '=*^PG)|| 

\r<i IIl2(7.v) 

EOC 

0 

4.0298e-02 

— 

8.2000e-01 

— 

1 

1.1316e-02 

1.83 

2.2827e-01 

1.84 

2 

1.4371e-03 

2.98 

2.8876e-02 

2.98 

3 

1.8037e-04 

2.99 

3.6208e-03 

3.00 

4 

2.2569e-05 

3.00 

4.5295e-04 

3.00 

5 

2.8219e-06 

3.00 

5.6631e-05 

3.00 


Table 5.2: 


Norm values and corresponding experimental orders of convergence in space- 
time for cGP(2)-MFEM(2) on the refinement levels as given in Tab. 


5.1 


domain with polynomial order r = 2 and p = 2, respectively, in the definition of the 
underlying finite element spaces. We note that the optimal order convergence in time 
and space is obtained for the primal and the flux variable. Thus, the estimates (4.30) 


and (4.29) might be suboptimal with respect to the time discretization; cf. Rem. 4.7 


The estimate (4.43) is nicely confirmed by the presented numerical results. Further, 
we note that the optimal rate of convergence is obtained for the spatial discretization 
of the flux field in the norm of V. In this point the family of Raviart-Thomas pairs 
of mixed finite elements is superior to the family of Brezzi-Douglas-Marini pairs of 
mixed finite elements (cf. [T7]) for that the optimal order of convergence of the flux 
variable can be obtained only in the norm of 


5.2. Distorted meshes 

In the second part of the numerical convergence studies we approximate the same 
analytic solution as before but we use spatial meshes with randomly distorted interior 
vertices. Precisely, each of the interior vertices is distorted by a randomly chosen 
vector. The magnitude of the distortion vector is chosen randomly up to a given 
factor of relative length to the corresponding edge length. The characteristic numbers 
of the refinement levels are summarized in Tab. |5.3| The resulting distorted meshes 
are illustrated in Fig. |5.2| for the refinement level 3. The temporal mesh is chosen in 
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Figure 5.1: Calculated errors and experimental orders of convergence in space-time for 
cGP(2)-MFEM(2). 


Level 

0% 

^max ^black 

5% 

^max ^black 

10% 

^max ^black 

25% 

^max ^black 

0 

1.4142 

— 

1.4142 

— 

1.4142 

— 

1.4142 

— 

1 

0.7071 

2.00 

0.7258 

1.95 

0.7504 

1.88 

0.8127 

1.74 

2 

0.3536 

2.00 

0.3661 

1.98 

0.3974 

1.89 

0.4694 

1.73 

3 

0.1768 

2.00 

0.1888 

1.94 

0.2006 

1.98 

0.2393 

1.96 

4 

0.0884 

2.00 

0.0946 

2.00 

0.1008 

1.99 

0.1195 

2.00 


Table 5.3: Distorted spatial mesh: /imax largest cell diameter and hbiack cell diameter 
reduction factor for 0%, 5%, 10% and 25% random vertex movement. 


the same way as in the first numerical experiment; cf. Tab. |5.l] 

We summarize the calculated errors and the corresponding experimental order of con¬ 
vergence (EOC) for the proposed space-time discretization on the distorted spatial 
meshes in Tab. |5.4]for the scalar-valued primal variable and in Tab. |5.5]for the vector¬ 


valued flux variable and further illustrate them in Fig. 5.3 The expected experimen¬ 
tal order of convergence in space and time, for the primal variable measured in the 
L^(/; L^(D))-norm and for the flux variable in the L^{I\ V)-norm, is largely confirmed 
even for the strongly perturbed meshes with a distortion factor of 25 %. This nicely 
demonstrates the robustness of the numerical scheme. We note that the space-time 
convergence studies on the distorted spatial meshes were done with exactly the same 
numerical solver settings as for the above-given studies on uniform meshes. 
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(a) 5% 


(b) 10% 



Figure 5.2: Distorted spatial meshes for 5% (a), 10% (b) 
movement for refinement level 3. 


and 25% (c) random vertex 


Level 

5% 

EOC 

10% 

EOC 

25% 

EOC 

0 

4.0298e-02 

— 

4.0298e-02 

— 

4.0298e-02 

— 

1 

1.1353e-02 

1.83 

1.1463e-02 

1.81 

1.2155e-02 

1.73 

2 

1.4381e-03 

2.98 

1.4783e-03 

2.95 

1.8690e-03 

2.70 

3 

1.8290e-04 

2.97 

1.9117e-04 

2.95 

2.5232e-04 

2.89 

4 

2.3024e-05 

2.99 

2.4402e-05 

2.97 

3.4463e-05 

2.87 


Table 5.4: Calculated errors and corresponding experimental ord er of convergence for 
11^2(/. 2,2 on distorted meshes given in Tab. 


5.3 


Level 

5% 

EOC 

10% 

EOC 

25% 

EOC 

0 

8.2000e-01 

— 

8.2000e-01 

— 

8.2000e-01 

— 

1 

2.2964e-01 

1.84 

2.3376e-01 

1.81 

2.6241e-01 

1.64 

2 

3.0172e-02 

2.93 

3.5136e-02 

2.73 

7.5321e-02 

1.80 

3 

3.9244e-03 

2.94 

4.7838e-03 

2.88 

1.0391e-02 

2.86 

4 

6.0990e-04 

2.69 

9.6250e-04 

2.31 

2.7408e-03 

1.92 


Table 5.5: Calculated errors and corresponding experimental orde r of convergence for 
||^cGP( 2 ) distorted meshes given in Tab. 


5.3 
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Figure 5.3: Calculated errors and corresponding experimental order of convergence on 
distorted meshes given in Tab. 


5.3 


6. Conclusions 

In this work a numerical analysis of a family of variational space approximation 
schemes that combine continuous finite elements in time with the MFEM in space 
was presented for a parabolic prototype model of flow in porous media. The existence 
and uniqueness of the temporally semidiscrete and the fully discrete approximations 
were proved. Error estimates with explicit rates of convergence, including an optimal 
order error estimate, in natural norms of the scheme were established. The error es¬ 
timates were illustrated and confirmed by numerical convergence studies. We believe 
that our analyses and techniques can be extended and applied to more sophisticated 
flow and transport processes in porous media or to incompressible viscous free flow. 
This will be our work for the future. 
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A. Supplementary proofs 


In the sequel we introduce a variational semidiscretization in time of the weak formu¬ 


lation of the second order problem (2.6)-(2.8), i.e. without rewriting Eq. (2.6) as a first 


order system of equations as it is done in Subsec. |2.3[ Then we prove the existence 
and uniqueness of solutions to the resulting semidiscrete variational problem. This 
results is used to establish the existence of the semidiscrete approximation in mixed 
form defined by the variational problem (2.14), (2.15) in Sec. Here, we present a 
different technique of proof than in |46) since one of the arguments that is used |461 
Lemma 6.1] does not hold in the applied form from our point of view. Thereby, we 
aim to keep our work self-contained. Further, we summarize the proof of Thm. |3.10[ 


A.l. Variational time discretization of the second order problem 

In the following we use the notation that is introduced in Subsec. [O] and [2731 respec¬ 


tively. Moreover we use the splitting (cf. Eq. (3.12)) 

Ur{t) = Uo +u^(t) with G Xq{Hq{Q)) . 
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Further, we put 


fo{t) = fit) - Auo 


2.1 


Under the 


for uo € H^i^l) such that by definition Auq € H ^(U), cf. Sec. 
additional regularity condition that 

DiA) = H^{n)nH^in) 

it even holds that /o € L^(/;L^(U)). The semidiscrete variational approximation of 
the system (2.6|-(2.8) is now defined by: Find S Xg(i/Q(r2)) such that 


[ {dtU°,Wr)dt+ f aiu°,Wr)dt= ( {fQ,Wr)dt (A.l) 

Jo Jo Jo 

for all Wr G (fl)). 

Firstly, we show the uniqueness of solutions to (A.l). In the sequel we denote by 
Fn,j = Fn,jit) for j = 0, ...,r the Lagrange basis functions in J„ = with 

respect to r + 1 quadrature points tn,i, I = 0, ...,r. Here, we choose the Gauss- 
Lobatto quadrature rule that is exact for polynomials of maximum degree 2r — 1. In 
particular, for the quadrature nodes in /„ it holds that tn,o = tn-i and tn^r = tn- 
Then, any function G X^iH^ifl)) and its time derivative admit the representation 


Urit) = UnFn,jit) , ^ U^lf'^jit) 

j=0 j=0 

for all t G In with coefficient functions G Fig (H) for j = 0,..., r. 


(A.2) 


Theorem A.l (Uniqueness of solutions to ( |A.1[ )). Let the assumptions of Subsec. 2.2 
about Lt,uo and f he satisfied. Then the solution G Aq (iFQ(H)) of the semidiscrete 
problem (A.l) is unique. 


Proof. Let 2 G '^o(Flo(U)) denote two solutions of the semidiscrete vari¬ 
ational problem (A.l). We put u^{t) := ^ — it° 2 - We choose the test function 

Wt '■= A~^dtu^ + for some fixed parameter /r > 0. By means of (A.2), it holds 

that Wr G (U)). For this choice of Wr it follows that 


I := 


.T 

(i9tu°(t),A ^dtUr + fidtUr) dt + {Au°it),A -l-/r9tu°) dt = 0. (A.3) 

Jo 


By the symmetry of a(-, •) we have that a(M°, dtUr) = j m°). Further, we have 

that {Ur,dtUr) = ysll'«?lli 2 (Qp Recalling (2.3)-(2.5l and noting that m°( 0) = 0 and 
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dtu° G it follows from (A.3) that 


0 = / = 


((9tu°(t),A '^dtu°)dt+ [ {dtu°{t),ndtu°) dt 
JO 

+ [ {Aur{t),A~^dtu°) dt + f {Au°{t),fidtu°)dt 
Jo Jo 

>c f + / ||9tW°||^2(n)dt 

Jo Jo 


1 d 


> c / \\dtu. 
Jo 


°llff-i(n) 


dt + ^ 


Il'9tw?||i2(a) dt 


+ 2 ll“T(^)lli2(n) + ■ 

This implies that = 0 and, consequently, that ^ = v 


The uniqueness of 


solutions to (A.l) is thus established. 


We remark that testing Eq. (A.l I with Vr = A~^dtu'^ or Vr = drU° would already be 
sufficient for proving the uniqueness result. Further, the symmetry of a(-, •) is essential 
in the previous proof. A generalization of the arguments to problems with nonsym- 
metric bilinearforms, for instance to convection-diffusion equations, still remains an 
open problem. 


The existence of a solution to the semidiscrete problem (A.l) follows from the unique¬ 


ness of the solutions. Using the eigenspaces of A, problem (A.l) can be reduced to 


a set of finite dimensional problems, for each of which obviously uniqueness implies 
existence. For this we recall the following result from Appendix D.6]. 

Lemma A.2. Let H be a separable Hilbert space, and suppose that S : H i-G H is a 
compact and symmetric operator. Then there exists a countable orthonormal basis of 
H consisting of eigenfunctions of S. 


Theorem A. 3 (Existence of solutions to (A.l)). Let the assumptions of Subsec. 2.2 


about fl,UQ and f as the be satisfied. Then the semidiscrete problem (A.l) admits a 
solution uf € Xq{Hq{LI)). 


Proof. The operator S := A~^ : L‘^{il,) L^{Tt) with A being defined in (2.1) is 
a bounded, linear compact operator mapping L^{Ll) into itself. By means of Lemma 
A.2 there exists a set of appropriately scaled eigenfunctions {wk}^i C L'^{Ll) with 
Wk G Hq{LI) such that {wk}^-^ is an orthogonal basis of Ffo(U) and an orthonormal 
basis of L'^{Vt). 
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In terms of these eigenfunctions {wk}^^i C Hq{Q) the solution of problem 
can be represented as 

r r oo 

Ur{x,t) , for t G 7„ . 

j—Q j—Ok—1 

(i) 

with coefficients d^), G K for fc = 1,..., oo and each j = 0, ..., r and n = 0, ..., iV. 
We choose test functions Vr € (fl)) being defined by 

Wkipii , for t G, 

0 , for t G I\In 

for * = 1,..., r, k = 1,..., oo and n — 1,... ,N. Then, for each k = 1,..., oo, we get 
the finite dimensional problem 



A.l) 




3=0 


3=0 


:=7fc 


= A.: 


Jo 


(A.4) 


= -bk,i 


for i = 1,... r and n = 1,... ,N. Due to the continuity of functions Ur G Ag (-ffQ(fl)) 
and the choice of the Gauss-Lobatto quadrature rule it holds that 


Tj(0) _ Tj{r) 
'^n — '^n-1 


or 


0'n,k ~ “n-l,/c ’ 


respectively. Therefore, we recast the finite dimensional problem (A.41 as 


{a.j + -fkPij) = h,i + {oiio + 7feAo) 


(A.5) 


for i = 1,... r and each k = 1,..., oo. For the finite dimensional problem (A.51 the 
uniqueness of the solution established in Thm. ED then implies the existence of a 
solution. ■ 


Similarly to Corollary |3.4[ the existence and uniqueness of the semidiscrete solution 
implies that an inf-sup stability condition in the underlying space-time framework is 
satisfied [Ml p. 85, Thm. 2.6]. 

Corollary A.4. Let the assumptions of Subsec. \2.^ about Lt,uo,D and f be satisfied. 
Letu^ G Xq{Hq{Q)) be the unique solution of the semidiscrete problem (A.l) according 
to Thm. m and\A.!^ Then, there exists a constant c > 0 such that 


dd^U^., V.p') 


inf 


> c> 0 
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with 


B{Ur,Vr)= / {dtUr + AUt , Vr) dt , 

Jo 


\Ut\\x — 


1/2 


(A.6) 


llfrlli^ = l|Wr||L2(/,ffi(n)) • 


Remark A.5. By the arguments of Thm. 6.2] the inf-sup stability condition 
implies an error estimate for the semidiscretization where the error is measured 

in the corresponding natural norm (A.6) of the scheme. 


A.2. Proof of Thm. [TTOl 

Proof. We let Wt := ItU — Ur, = J^q — q^. with the interpolation operators Ir 


and Jr of Subsec. 3.2 Using the inf-sup stability condition along with problem (2.10) 


(2.11) and the non-exact semidiscrete problem (3.39), (3.40), applying the inequality 


of Cauchy-Schwarz and the continuity (2.13) of ar{-, •) we conclude that 

a\\{Wr,Vr}\\w\\{^r'lfr}\\v < ar{{Wr ,V r}, {ffr r}) 

= a^({u°, q}, {ipr, '0r}) - ar{{u°r, q^}, {ifr, if r}) 

- ar{{u° - Uu°, q - J^q}, {(Pr, ifr}) 

= [ {f{t)-Arfit),Pr)dt-ari{u° -IrU°,q-Jrq},{fr,lfr}) 

Jo 

< 11/ - Arf\\L^(I-,W)\\Wr,lfr}\\v + c||{u° - IrU°,q- J rqjWwWWr , ll’r}\\v 

( ^ 

— ( ^ 111/ ~ 11^/11^2(7.ly) -I-c||u° — 

\ 

\ 1/2 

+ c||5t(u° -/^U°)||^2(7_^.iy) +c||q- J^q||^2(7^.V)} j ’ l!{</’r, ^/’rlllv ■ 


By means of the approximation properties (3.16) to ( 3.18[ ) we then get that 

/ N 

\\{Wr,Vr}\\w < ci E^n’'{l|9r^||i2(7„;n.) 

\ n^l 

+ r^n\\dl+^q\\ 

Combining this estimate with the triangle inequality yields the assertion of Thm. |3.10| 


1/2 
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B. Summary of notation 

Function spaces and norms 


I, Q 

time and space domain, I = (0, T] 

L\n), 

standard Sobolev spaces 

H-\Q) 

dual space of Hq{Q.) 

1! ■ li> II ■ llpi ('i■) 

Norms in inner product in L^{Q) 

W, V 

W = L2(n), V = H{dw;n) 

II ' llv 

ll•lk = (||•f + ||v.(.)f)'/" 

C{7;X), L\I-X), H\I-X) 

Bogner spaces with values in the Banach space X 

D, uo, f 

data of the problem 

u, q 

solution and flux of the continuous problem 

a{u, v) 

bilinear form a(u, v) = (Du, v) 


Time discretization 



trial spaces of semidiscretization in time 

XS(W) 

trial space of semidiscretization with Wt(0) = 0 

y^(v) 

test spaces of semidiscretization in time 

Pr-(J; X) 

p.(J; X) = {p{t) = a 1 e e X} 

W, V 

W = XS(W) X X^V), V = Y'-^(W) X Y''-^(V) 

-fn, Tn 

subinterval 7„ = t„ = t„ — t„-i 

^n,0 

^n,0 — in — 1 

in,l ? • • • ) i'n,r 

r-point Gaussian quadrature nodes 

‘Pn,j{t) 

Lagrange polynomial on I„ w.r.t t„,o, • • •, tn,r 


Lagrange polynomial on /„ w.r.t t„,i,..., tn,r 

Uq- 

Vn,j(t) 

semidiscrete solution 

Qt 

<?T|I„(i) = Ei=0 Qn ‘PnoW 
semidiscrete flux 

! 

I = [0,1] reference interval 

... ^tr 

Gaussian quadrature nodes on I 

a;i,..., c2?r 

Gaussian quadrature weights on I 


transformed Lagrange polynomial on I 

&ij, Pij 

&ij = (hi, Pij = hiSij with Kronecker symbol 5ij; 
from transformation of time integrals to I and 
application of quadrature on I 

Hr 

temporal Lagrange interpolant 
w.r.t. to tnfi, ... 5 tn,r 

Ir 

temporal interpolation operator for variable u 

Jr 

temporal interpolation operator for variable q 

ar{{; ■},{;■}) 

space-time bilinear form 

u°{t), u°{t) 

U°(t) = U{t) — Uo, Ur{t) = Urit) — Uo 

z, p 

solution of dual problem 

h 

temporal interpolation operator for 2 

Jo 

temporal interpolation operator for p 
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Space discretization and error analysis 


Wh, Vh 

inf-sup stable pair of finite element spaces 
Raviart-Thomas (-Nedelec) elements 

X^{Wh), X^{Vh) 

trial spaces of space-time discretization 

y^-\Wh), y'-^Vh) 

test spaces of space-time discretization 

UtjH 

fully discrete solution 

Qr,h 

4'T.hW|i„ = Ej=0 Qi,hV>n,j{t) 
fully discrete flux 

Ph, Ph 

projection onto Wh and Vh 


projection onto Vh’- 

{V • (Iljin — v), Wh) = 0 for all Wh € Wh 

Eu{t) 

Eu(i) — Uxi^t) Ur,h(i) 

Eu{t) = Ei,n‘Pn,j{t) 

error between semidiscrete and fully discrete 

approximation of scalar variable 

E^{t) 

E<j{t) ~ Ej=0 E'q,'n‘fin,j{t) 

error between semidiscrete and fully discrete 

approximation of flux variable 

K,n 

Eu^ri = Eu{tn,i)- error in node t„,i, i = 0,..., r 


= Eq{tn,i)- crror in node t„,i, i = 0,..., r 
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